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North Atlantic climate during glacial times was characterized by large-amplitude 
switchings, the Dansgaard-Oeschger (DO) events, with an apparent tendency 
to recur preferably in multiples of about 1470 years. Recent work interpreted 
these intervals as resulting from a subharmonic response of a highly nonlin- 
ear system to quasi-periodic solar forcing plus noise. This hypothesis was chal- 
lenged as inconsistent with the observed variability in the phase relation be- 
tween proxies of solar activity and Greenland climate. Here we reject the claim 
of inconsistency by showing that this phase variability is a robust, generic 
feature of the nonlinear dynamics of DO events, as described by a model. 
This variability is expected from the fact that the events are threshold cross- 
ing events, resulting from a cooperative process between the periodic forc- 
ing and the noise. This process produces a fluctuating phase relation with 
the periodic forcing, consistent with proxies of solar activity and Greenland 
climate. 
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1. Introduction 

Climate archives from the North Atlantic region show repeated shifts in glacial climate, 
the Dansgaard-Oeschger (DO) events [Grootes et al., 1993; Andersen et al., 2006]. During 
Marine Isotope Stages 2 and 3 the intervals between the events exhibit a tendency to co- 
incide approximately with multiples of 1470 years [Rahmstorf , 2003] , as depicted in figure 
1. The statistical significance of this pattern and the responsible mechanism, however, 
is still a matter of debate [Ditlevsen et al., 2005]. Several hypotheses were proposed to 
explain the timing of DO events. One of these relates the events to two century-scale 
solar cycles with periods close to 1470/7 (=210) and 1470/17 («87) years [Braun et al, 
2005], the so-called De Vries/Suess [Wagner er al., 2001] and Gleissberg [Peristykh and 
Damon, 2003] cycles. Support for a leading solar role comes from deep-sea sediments, 
which indicate that during the Holocene century-scale solar variability was a main driver 
of multi-centennial scale climate changes in the North Atlantic region [Bond et al., 2001]. 

Recently the phase relation between solar variability (deduced from 10 Be) and 14 DO 
events was analyzed [Muscheler and Beer, 2006]. A relation far from fixed was found and 
was interpreted as being in contradiction to Braun et al.'s hypothesis [Braun et al., 2005]. 
While in linear systems a constant phase relation between the forcing and the response 
is expected, such a relation does not necessarily exist in non-linear systems. But climate 
records and ocean-atmosphere models [Ganopolski and Rahmstorf, 2001], which are not 
yet suitable for statistical analyses on DO events because of their large computational 
cost, suggest that the events represent switches between two climate states, consequently 
implying an intrinsically non-linear dynamical scenario. Thus, to interpret the reported 
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lack of a fixed phase relation between the DO events and solar proxies, it is crucial to 
analyze their phase relation in simple models. 

2. A Simple Model of DO Events 

Here we investigate this phase relationship in a very simple model of DO events. A 
comprehensive description of this model has already been published before, including a 
detailed discussion of its geophysical motivation and its applicability, as well as a com- 
parison with a much more detailed ocean-atmosphere model [Braun et al., 2007]. In the 
simple model, which was derived from the dynamics of the events in that ocean-atmosphere 
model, DO events represent repeated switches between two possible states of operation 
of a bistable, excitable system with a threshold (Figure 2). These states correspond to 
cold and warm periods of the North Atlantic region during DO cycles. The switches are 
assumed to occur each time the forcing function (f), which mimics the solar role in driving 
DO events, crosses the threshold (T). Transitions between the two states are accompanied 
by an overshooting of the threshold, after which the system relaxes back to its respective 
equilibrium following a millennial-scale relaxation. 

The rules for the transitions between both states are illustrated in figure 2. It is assumed 
that the threshold function T is positive in the interstadial ("warm") state and negative 
in the stadial ("cold") state. A switch from the stadial state to the interstadial one is 
triggered when the forcing / is smaller than the threshold function, i.e. when f(t) < T(t). 
The opposite switch occurs when f(t) > T(t). During the switches a discontinuity in 
the threshold function is assumed, i.e. T overshoots and takes a non-equilibrium value 
(A during the shift into the stadial state, Ai during the opposite shift). Afterwards, T 
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approaches its new equilibrium value (B in the stadial state, Bi in the interstadial state) 
following a millennial scale relaxation process: 

dT T — B< 



's 



(1) 



dt r s 

Here, r and T\ represent the relaxation time in the stadial state (s=0) and in the inter- 
stadial state (s=l), respectively. 

Both the overshooting relaxation assumption and the transition rules in our simple 
model are a first order approximation of the dynamics of DO events in a coupled ocean- 
atmosphere model [Ganopolski and Rahmstorf, 2001]. In that model the events also 
represent threshold-like switches in a system with two possible states of operation (cor- 
responding to two fundamentally different modes of deep water formation in the North 
Atlantic) and with an overshooting in the stability of the system during these shifts 
[Ganopolski and Rahmstorf, 2001; Braun et al., 2007]. Analogous to the simple model, 
switches from the stadial mode into the interstadial one are triggered by sufficiently large 
negative forcing anomalies (i.e. by a reduction in the surface freshwater flux to the North 
Atlantic that exceeds a certain threshold value), whereas the opposite shifts are trig- 
gered by sufficiently large positive forcing anomalies (i.e. by an increase in the freshwater 
flux that exceeds a certain threshold value). It has further been demonstrated that the 
simple model is able to reproduce the timing of DO events as simulated with the ocean- 
atmosphere model, as well as the occurrence of non-linear resonance phenomena such as 
stochastic resonance and ghost resonance, which were shown to be properties exhibited 
by that model [Ganopolski and Rahmstorf, 2002; Braun et al., 2005; Braun et al., 2007]. 
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An obvious advantage of the conceptual model compared with the ocean-atmosphere 
model is its low computational cost, which allows for extensive statistical analyses on the 
timing of DO events. All model-parameter values chosen here are the same as in two 
earlier publications (A = -27 mSv, Ai = 27 mSv, B = -9.7 mSv, Bi = 11.2 mSv, 
r = 1200 years, T\ = 800 years; 1 mSv = 1 milli-Sverdrup = 10 3 m 3 /s) [Braun et al., 
2005; Braun et al, 2007]. 

3. Testing Fixed Phase Relationships 

To test the assumption of a fixed phase relationship between solar-forced DO events 
and solar variability, we here drive our model by a simple input consisting of noise a ■ n(t) 
and of two sinusoidal cycles with equal amplitudes A: 

/(f) = -A ■ (cos® + cos®) + a ■ nit). (2) 

J-l J-2 

a is the standard deviation of the noise and n(t) the standard unit variance white noise, 
with a cutoff frequency of 1/50 years (figure 3). Following Braun et al. [2007] the cutoff 
is used to account for the fact that the model shows an unrealistically large sensitivity to 
decadal-scale or faster forcing. In analogy to Braun et al. [2005] the periods of the two 
cycles are chosen to be T x = 1470/7 (=210) years and T 2 = 1470/17 («86.5) years, i.e. 
close to the leading spectral components of the solar De Vries and Gleissberg cycles. 

In the simulations shown in Figure 4 we use three different signal-to- noise ratios (SNR): 
A = 8 mSv and a = 5.5 mSv (SNR « 2.1), A = 5 mSv and a = 8 mSv (SNR « 0.4), A 
= 3 mSv and a = 9 mSv (SNR rs 0.1). For all of these, the waiting time distribution of 
the simulated events is centered around a value of 1470 years, with several peaks of only 
decadal-scale width (figure 4). The relative position of these peaks is well understood 
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in the context of the ghost resonance theory [Chialvo et al., 2002; Chialvo, 2003; Calvo 
and Chialvo, 2006; Braun et al., 2007]. The peaks result from constructive interference 
between the two sinusoidal forcing cycles which produces particularly large magnitude 
variations in the bi-sinusoidal forcing and - when noise is added - leads to favored tran- 
sitions at the corresponding waiting times. Depending on the relative amplitude values 
of the noise and the periodic forcing this synchronization is more or less efficient, as is 
seen for the different signal-to-noise ratios in Figure 4. Even for the lowest ratio, however, 
the synchronization is still notable. The waiting time distributions shown in figure 4 are 
thus almost symmetrically centered around a preferred value of 1470 years because the 
sinusoidal cycles enter in phase every 1470 years, creating forcing peaks of particularly 
large magnitude. This 1470-year repeated coincidence of the bi-sinusoidal forcing, how- 
ever, does not show up as a corresponding forcing frequency, since no sinusoidal cycle with 
that period is present. Thus, when linear spectral analysis is performed on the forcing, 
only the two century-scale sinusoidal cycles are detected as outstanding components. 

Despite the robustness of the synchronization effect, none of the two sinusoidal cycles 
in our forcing shows a fixed phase relationship with all of the simulated DO events, due to 
the presence of noise and the existence of a threshold. In our model, a fixed phase relation 
can only be present in the low noise limit (i.e. either for a — > [with a supra-threshold 
bi-sinusoidal forcing] or for the lowest noise level that still enables repeated threshold 
crossings [with a sub-threshold bi-sinusoidal forcing], thus corresponding to DO events 
with extremely long waiting times), compare third column in figure 4. Even for the largest 
of the three signal-to-noise ratios in our simulations, the events thus only show a tendency 
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to cluster around a preferred phase of the two sinusoidal cycles. Outliers, however, can still 
occur in almost opposite phase, at least once over a sufficiently large number of events in 
the simulation. For the highest signal-to-noise ratio, for example, there is still a probability 
of about 35 percent to find at least one out of 14 events in opposite phase (i.e. outside of 
the interval [— vr/2, +7r/2]). And for the other two signal-to-noise ratios, the corresponding 
probability is even much higher (i.e. 92 percent and 99.5 percent, respectively). Since a 
fixed phase relationship between the simulated events and the forcing cycles does not even 
exist in our very simple model system, it appears unrealistic to us to assume the existence 
of such a relationship in the climate system. Thus, the reported lack of a fixed phase 
relationship between DO events and solar variability (deduced from 10 Be) would also be 
expected with the proposed ghost resonance solar forcing. We note that superimposed 
epoch analyses of 14 simulated events can produce forcing-response relations similar to 
the one reported by Muscheler and Beer [2006]: The onset of the superimposed events 
(at t — in the fourth column in figure 4) typically coincides with a minimum in the 
averaged bi-sinusoidal forcing which, however, is not more pronounced than other minima 
and is highly damped as compared with the unaveraged forcing. A considerable statistical 
spread exists in the magnitude of this damping because the small number (14) of events 
is not yet sufficient to infer reliable information concerning the average phase fluctuation 
between the input and the output. 

This lack of phase correlation between forcing and response is explained by the threshold 
character of DO events: The simulated events are triggered when the total forcing (the 
sum of the two sinusoidal cycles and the noise) crosses the threshold function. Some 
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of the threshold-crossings are in the first place caused by constructive interference of 
both cycles. These events coincide with near-minima of the two forcing cycles. Other 
threshold-crossings are, however, in the first place caused by constructive interference 
of just one cycle and noise. These events thus coincide with a near-minimum of only 
that cycle (compare figure 3), whereas a fixed phase-relation with the second cycle does 
not necessarily exist. And at least for low signal-to-noise ratios, some of the threshold 
crossings are in the first place caused by the noise alone. These events thus do not show 
a fixed phase relation with any of the two forcing cycles. 

The inherently nonlinear noisy synchronization mechanism exhibited by our model is not 
unique to DO events. In fact, it has originally been proposed to explain the perception of 
the pitch of complex sounds [Chialvo et al., 2002; Chialvo, 2003] and, as a general concept, 
has already been used to describe theoretically and experimentally similar dynamics in 
other excitable [Calvo and Chialvo, 2006] or multi-stable systems with thresholds, e.g. in 
lasers [Chialvo et al., 2002; Chialvo, 2003; Buldu et al., 2003]. Because of the fact that 
the leading output frequency is absent in the input, this type of resonance is called ghost 
stochastic resonance. 

4. Conclusions 

We here used a simple model of DO events, driven by a bi-sinusoidal forcing plus noise, 
to show that a fixed phase relation between the forcing cycles and all simulated events 
does not exist, apart from the unrealistic low noise limit. As argued above, in this model 
the fluctuations in the phases between the forcing and the response are related to the 
process giving rise to the transition itself. Each event is generated by a threshold crossing 
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resulting from a cooperative process between the two periodic driving forces (i.e., the 
centennial-scale input cycles) and the stochastic fluctuations. In this nonlinear scenario, 
as we showed explicitly in our simulations, millennial-scale events with fixed input-output 
phase relations are impossible for any nonzero noise amplitude. 

While one could disagree on the interpretation and the statistical significance of the pat- 
tern described in Figure 1, our results show that the reported lack of a fixed phase relation- 
ship between 14 DO events and solar proxies is consistent with the suggested solar role in 
synchronizing DO events. At the same time our results have further implications for a sec- 
ond so-far unexplained oscillation during Pleistocene climate, i.e. the glacial-interglacial 
cycles, which also show strong indications for the existence of threshold-like dynamics dur- 
ing glacial terminations [Paillard, 1998; Huybers and Wunsch, 2005]: Since the existence 
of a causal relation between threshold-crossing events and their quasi-periodic forcing 
does not necessarily imply the existence of a clear phase relation over all events, the lack 
of such a relation between glacial terminations on one hand and the orbital eccentricity 
and precession cycles on the other hand is not sufficient to reject a leading role of these 
cycles during terminations, in contrast to the interpretation proposed by Huybers and 
Wunsch [2005]. More insight in the cause of Pleistocene climate cycles might thus be 
gained from more adequate statistical approaches, based e.g. on Monte-Carlo simulations 
[Ditlevsen et al., 2007] with simple models [Paillard, 1998; Braun et al., 2007] that mimic 
the nonlinear dynamics which seems to be relevant during these oscillations. 



DRAFT 



March 31, 2008, 11:56am 



DRAFT 



BRAUN, DITLEVSEN AND CHIALVO: CLIMATE SHIFTS & LINEAR PHASE RELATIONS X - 11 
Acknowledgments. H.B. thanks A. Ganopolski for instructive discussions on the dy- 
namics of DO events and A. Svensson for his hospitality at the Niels Bohr Institute. H. 
B. was funded by Deutsche Forschungsgemeinschaft (project MA 821/33). 

References 

Andersen, K., A. Svensson, S. J. Johnsen, S. O. Rasmussen, M. Bigler, R. Rothlisberger, 

U. Ruth, M.-L. Siggaard-Andersen, J. P. Steffensen, D. Dahl- Jensen, B. M. Vinther, 

and H. B. Clausen (2006), The Greenland Ice Core Chronology 2005, 15-42 ka. Part 1: 

Constructing the time scale. Quat. Sci. Rev., 25, 3246-3257. 
Bond, C; B. Kromer, J. Beer, R. Muscheler, M. N. Evans, W. Showers, S. Hoffmann, 

R. Lotti-Bond, I. Hajdas, and G. Bonani (2001), Persistent Solar Influence on North 

Atlantic Climate During the Holocene, Science, 294, 2130-2136. 
Braun, H., M. Christl, S. Rahmstorf, A. Ganopolski, A. Mangini, C. Kubatzki, K. Roth, 

and B. Kromer (2005), Possible solar origin of the 1,470-year glacial climate cycle 

demonstrated in a coupled model, Nature, 438, 208-211. 
Braun, H., A. Ganopolski, M. Christl, and D. R. Chialvo (2007), A simple conceptual 

model of abrupt glacial climate events, Nonlin. Processes Geophys., 14, 709-721. 
Buldu, J. M., D. R. Chialvo, C. R. Mirasso, M. C. Torrent, and J. Garcia-Ojalvo (2003), 

Ghost resonance in a semiconductor laser with optical feedback, Europhysics Letters, 

64, 178-184. 

Calvo, O. and D.R. Chialvo (2006), Ghost stochastic resonance on an electronic circuit. 
International Journal of Bifurcation and Chaos, 16, 731-735. 



DRAFT 



March 31, 2008, 11:56am 



DRAFT 



X - 12 BRAUN, DITLEVSEN AND CHIALVO: CLIMATE SHIFTS & LINEAR PHASE RELATIONS 

Chialvo, D. R., O. Calvo, D. L. Gonzalez, O. Piro, O, and G. V. Savino (2002), Subhar- 
monic stochastic synchronization and resonance in neuronal systems, Physical Review 
E, 65, 050902. 

Chialvo, D. R. (2003), How we hear what is not there: A neuronal mechanism for the 

missing fundamental illusion, Chaos, 13, 1226-1230. 
Ditlevsen, P. D., M. S. Kristensen, and K. K. Andersen (2005), The recurrence time of 

Dansgaard-Oeschger events and limits on the possible periodic component, Journal of 

Climate, 18, 2594-2603. 
Ditlevsen, P. D., K. K. Andersen, and A. Svensson (2007), The DO-climate events are 

probably noise induced: statistical investigation of the claimed 1470 years cycle, Clim. 

Past, 3, 129-134. 

Ganopolski, A. and S. Rahmstorf (2001), Simulation of rapid glacial climate changes in a 
coupled climate model, Nature, 409, 153-158. 

Ganopolski, A. and S. Rahmstorf (2002), Abrupt glacial climate changes due to stochastic 
resonance, Phys. Rev. Lett., 88, 038501. 

Grootes, P. M., M. Stuiver, M., J. W. C. White, S. Johnsen, and J. Jouzel (1993), Compar- 
ison of oxygen isotope records from the GISP2 and GRIP Greenland ice cores, Nature, 
366, 552-554. 

Huybers, P. and C. Wunsch (2005), Obliquity pacing of the late Pleistocene glacial ter- 
minations, Nature, 434, 491-494. 

Muscheler, R. and J. Beer, Solar forced Dansgaard/Oeschger events (2006)? Geophys. 
Res. Lett, 33, L20706. 



DRAFT 



March 31, 2008, 11:56am 



DRAFT 



BRAUN, DITLEVSEN AND CHIALVO: CLIMATE SHIFTS & LINEAR PHASE RELATIONS X - 13 

Paillard, D. (1998), The timing of Pleistocene glaciations from a simple multiple-state 

climate model, Nature, 378, 378-381. 
Peristykh, A. N. and P. E. Damon (2003), Persistence of the Gleissberg 88-yr solar cycle 

over the last 12,000 years: Evidence from cosmogenic isotopes. J. Geophys. Res., 108; 

1003. 

Rahmstorf, S. (2003), Timing of abrupt climate change: a precise clock, Geophys. Res. 
Lett, 30, 1510-1514. 

Wagner, G., J. Beer, J. Masarik, R. Muscheler, P. W. Kubik, W. Mende, C. Laj, G. M. 
Raisbeck, and F. Yiou (2001), Presence of the solar de Vries cycle (^205 years) during 
the last ice age, Geophys. Res. Lett., 28, 303-306 



DRAFT 



March 31, 2008, 11:56am 



DRAFT 



X - 14 



BRAUN, DITLEVSEN AND CHIALVO: CLIMATE SHIFTS & LINEAR PHASE RELATIONS 




O cold L^l — ! — 1 — ! — , — , — !_J — i — 1 — , — , — 1 — — 1 — , — , — , — i — , 

40000 30000 20000 10000 

Time (years before present) 



Figure 1. DO events as seen in two ice cores from Greenland. Top: NGRIP. Bottom: 
GISP2. Labels mark the events 1-10, the Aller0d (A) and the end of the Younger Dryas 
(0). Dashed lines are spaced by 1470 years. 
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Time 

Figure 2. Dynamics of DO events in the model. Top: Forcing f (grey) and threshold 
function T (black). Bottom: Model state s (s=0: "cold" state, s — 1: "warm" state). 
A switch from the cold to the warm state is triggered when / < T, which happens in 
this example at time t" . During the transition, interpreted as the start of a DO event, 
T overshoots and relaxes back towards its new equilibrium Bi (B\ > 0) following a 
millennial time scale. The events are terminated by a switch back to the cold state, which 
is triggered when / > T (at time t' in the figure). Again, T overshoots and approaches 
its new equilibrium B (B < 0). 
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Figure 3. Forcing and response. The forcing consists of two added sinusoidal cycles 
with equal amplitudes (A = 8 mSv) and white noise (a = 5.5 mSv). (A) Total forcing 
f (grey) and threshold function T (black). (B) Forcing components (grey); from top to 
bottom: 210-year cycle, 86.5-year cycle, noise. Dashed lines indicate the onset of the 
simulated DO events. Despite the tendency of the three events to recur approximately 
every 1470 years, only the first two events coincide with minima of the 210-year cycle. 
The third event, in contrast, occurs closely after a maximum of that cycle. 
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Figure 4. Waiting time distribution of the simulated events and their phase relation 
with the forcing. The amplitude of the two sinusoidal forcing cycles and the standard 
deviation of the noise are: A = 8 mSv, a = 5.5 mSv (top); A = 5 mSv, a = 8 mSv 
(middle); A = 3 mSv, a = 9 mSv (bottom). First column: normalized distribution of 
the spacing between successive DO events in the simulation. Second column: probability 
distribution of the phase relation between the onset of the events and the 210-year cycle 
(zero corresponds to a start of the events at the minimum of that cycle, ±7r to a start at 
the maximum). Third column: Rayleigh's R value [Ditlevsen et al., 2007], a measure for 
the phase correlation between the forcing cycles and the simulated events, as a function 
of the noise level a (R = J2n=i cos + ^ • sin ^r 22 - 1 , where t n denotes the timing of the 
events, T\ = 210 years and iV — > oo). R is 1 if and only if a fixed phase relation exists 
between the 210-year forcing cycle and the events. Fourth column: Superposition of the 
model response (i.e. of the state variable s) and of the bi-sinusoidal forcing over a series 
of 14 simulated events, aligned by the onset of the events following Muscheler and Beer 
[2006]. The unaveraged bi-sinusoidal forcing is normalized with maximum and minimum 
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